##  function to estimate probit using MI variables
glm.MI <- function(d,outcome,add.pred.MI,add.pred,sims,wald)	{

	wrap <- function(x)	{ paste0("$(", formatC(x,format = "f",3), ")$") }

	out.coef	<- vector(mode = "list", length = 10)
	out.var		<- vector(mode = "list", length = 10)
	out.sims	<- vector(mode = "list", length = 10)
	out.n		<- vector(mode = "list", length = 10)

	out.converged <- rep(NA,10)


	v <- 1
	for(v in 1:10)	{

		d.temp <- d

		##  expand predictors that are multiply imputed
		add.pred.MI.expanded <- ""
		m <- 1
		for(m in 1:length(add.pred.MI))	{
			add.pred.MI.expanded <- paste0(add.pred.MI.expanded, " + ", gsub("@",v,add.pred.MI[m]))
										}

		f <- as.formula	(
		paste0(	outcome, " ~ ", add.pred.MI.expanded,
				" + FirstSessionYear + Denunciation + Euthanasia + FinalPhase +",
				"ExterminationEinsatzgruppen + ExterminationCamps + ExterminationOther +",
				"Judicial + NSDetainment + NSOther + WarCrimes",
				add.pred)
						)


		##  fit model twice to find out whether any obs/covars
		##  need to be dropped because of missingness
		glm.out <- glm(f, data = d.temp, family = binomial(link = "probit"))
		if(!is.null(na.action(glm.out)))	{ d.temp <- d.temp[-na.action(glm.out),] }

		glm.out <- glm(f, data = d.temp, family = binomial(link = "probit"))
		out.coef[[v]]	<- coef(glm.out)[!is.na(coef(glm.out))]
		out.var[[v]]	<- vcovCL(glm.out, fix = TRUE, cluster = cbind(d.temp$Court,d.temp$FirstSessionYear))
		##  clustered at level of assignment

		##  consistency check
		stopifnot(length(out.coef[[v]]) == nrow(out.var[[v]]))

		##  record n
		out.n[[v]] <- nrow(d.temp)

		##  MLE convergence check
		out.converged[v] <- glm.out$boundary == FALSE & glm.out$converged == TRUE


		##  Monte Carlo simulations of first difference
		if(!is.null(sims))	{

			m0.hi <- m0.lo <- model.matrix(glm.out)[,!is.na(coef(glm.out))]

			## consistency check
			stopifnot(colnames(m0.hi) == names(out.coef[[v]]))

			if(!is.logical(m0.hi[,(sims+1)]))	{
				m0.hi[,(sims+1)] <- quantile(m0.hi[,(sims+1)], probs = .75)
				m0.lo[,(sims+1)] <- quantile(m0.lo[,(sims+1)], probs = .25)
												}

			if(is.logical(m0.hi[,(sims+1)]))	{
				m0.hi[,(sims+1)] <- 1
				m0.lo[,(sims+1)] <- 0
												}

			out.hi <- matrix(NA,nrow(m0.hi),20000)
			out.lo <- matrix(NA,nrow(m0.hi),20000)

			M <- 1
			for(M in 1:20000)		{
				gamma.tilde	<- mvrnorm(n = 1, mu = out.coef[[v]], Sigma = out.var[[v]])
				out.hi[,M]	<- pnorm(m0.hi %*% gamma.tilde)
				out.lo[,M]	<- pnorm(m0.lo %*% gamma.tilde)
									}

			out.sims[[v]] <- apply(out.hi,2,mean) - apply(out.lo,2,mean)

							}

						}

	##  consistency checks
	stopifnot(sum(is.na(out.coef)) == 0)
	stopifnot(sum(is.na(out.var)) == 0)
	stopifnot(length(unique(sapply(out.coef,length))) == 1)


	##  summarize simulation results
	if(!is.null(sims))	{
		out.sims2 <- list(mean(unlist(out.sims)),sd(unlist(out.sims)))
						} else out.sims2 <- NA


	##  compute correct MI coefs and standard errors
	temp <- data.frame(matrix(NA,length(out.coef[[1]]),2))
	rownames(temp) <- names(out.coef[[1]])
	colnames(temp) <- c("coef","se")

	m <- 1
	for(m in 1:nrow(temp))	{

		out <- 
		pool.scalar(Q = unlist(lapply(out.coef, function(x)x[m])),
					U = unlist(lapply(out.var, function(x)x[m,m])),
					n = min(unlist(out.n)))

		temp[m,1:2] <- c(out$qbar, sqrt(out$t))

							}


	##  format results
	temp$z <- abs(temp[,1] / temp[,2])

	asterisks <- rep("", nrow(temp))
	asterisks[temp$z >= qnorm(0.975)]  <- "^{*}"
	asterisks[temp$z >= qnorm(0.995)]  <- "^{**}"
	asterisks[temp$z >= qnorm(0.9995)] <- "^{***}"

	temp[,1] <- paste0("$", formatC(round(temp[,1],3), format = "f", 3), asterisks, "$")
	temp[,2] <- wrap(round(temp[,2],3))


	##  Wald test
	Wald.out <- rep(NA,length(wald))

	if(!is.null(wald))	{

		##  judges
		Q <- matrix(0,1,length(out.coef[[1]]))
		Q[1,wald[[1]][1]+1] <- 1
		Q[1,wald[[1]][2]+1] <- -1
		Wald.out[1] <- unlist(MIwaldtest(qhat = out.coef, u = out.var, Cdes = Q, rdes = 0)$stat[4])

		##  prosecutors
		Q <- matrix(0,1,length(out.coef[[1]]))
		Q[1,wald[[2]][1]+1] <- 1
		Q[1,wald[[2]][2]+1] <- -1
		Wald.out[2] <- unlist(MIwaldtest(qhat = out.coef, u = out.var, Cdes = Q, rdes = 0)$stat[4])

						}

	return(list(temp,formatC(round(Wald.out,3), format = "f", 3), unlist(out.sims2), unlist(out.n)))

}







#

